#include <iostream>
#include <cmath>
#include <vector>
#include "DifferenceTable.h"

class Interpolation
{
    public:
        int PointNum;
        DifferenceTable *diffTable;
        std::vector<double> inputPoint;
        std::vector<double> inputValue;
    public:
        Interpolation(int n, std::vector<double> &x, std::vector<double> &y): 
            PointNum(n), inputPoint(x), inputValue(y), diffTable(new DifferenceTable[2*n]){};
        ~Interpolation(){delete [] diffTable;};
        virtual void get_DiffTable() = 0;
        virtual double getInterpolationValue(double x) = 0;
};

class NewtonFormula:public Interpolation
{
    public:
        NewtonFormula(int n, std::vector<double> &x, std::vector<double> &y): 
            Interpolation(n, x, y){};
        ~NewtonFormula(){};
        void get_DiffTable();
        double getInterpolationValue(double x);
};

class HermiteFormula:public Interpolation
{
    public:
        std::vector<double> inputDerivative;
    public:
        HermiteFormula(int n, std::vector<double> &x, std::vector<double> &y, std::vector<double> &dy): 
            Interpolation(n, x, y), inputDerivative(dy){};
        ~HermiteFormula(){};
        void get_DiffTable();
        double getInterpolationValue(double x);
};



void NewtonFormula::get_DiffTable(){
    DifferenceTable f[PointNum][PointNum];
    for(int i = 0; i < PointNum; i++){
        f[i][0] = DifferenceTable(i, i, inputValue[i]);
        if (i == 0) continue;
        for (int j = 1; j <= i; j++) {
            f[i][j] = DifferenceTable(f[i-1][j-1].get_head(), f[i][j-1].get_tail(), 
            (f[i][j-1].get_value() - f[i-1][j-1].get_value())/(inputPoint[f[i][j-1].get_tail()] - inputPoint[f[i-1][j-1].get_head()]));
        }
    }
    for(int i = 0; i < PointNum; i++){
        diffTable[i] = f[i][i];
    }
}

double NewtonFormula::getInterpolationValue(double x){
    NewtonFormula::get_DiffTable();
    double p = diffTable[0].get_value();
    double prod = 1.0;
    for (int i = 1; i < PointNum; i++) {
        prod *= (x - inputPoint[i-1]);
        p += diffTable[i].get_value() * prod;
    }
    return p;
}


void HermiteFormula::get_DiffTable(){
    DifferenceTable f[2*PointNum][2*PointNum];
    for (int i = 0; i < PointNum; i++) {
        f[2*i][0] = DifferenceTable(i, i, inputValue[i]);
        f[2*i+1][0] = DifferenceTable(i, i, inputValue[i]);
        f[2*i+1][1] = DifferenceTable(i, i, inputDerivative[i]);
        if (i == 0) continue;
        for (int j = 1; j <= 2*i; j++) {
            f[2*i][j] = DifferenceTable(f[2*i-1][j-1].get_head(), f[2*i][j-1].get_tail(), 
            (f[2*i][j-1].get_value() - f[2*i-1][j-1].get_value())/(inputPoint[f[2*i][j-1].get_tail()] - inputPoint[f[2*i-1][j-1].get_head()]));
            f[2*i+1][j+1] = DifferenceTable(f[2*i][j].get_head(), f[2*i+1][j].get_tail(),
            (f[2*i+1][j].get_value() - f[2*i][j].get_value())/(inputPoint[f[2*i+1][j].get_tail()] - inputPoint[f[2*i][j].get_head()]));
        }
    }
    for(int i = 0; i < 2*PointNum; i++){
        diffTable[i] = f[i][i];
    }
}

double HermiteFormula::getInterpolationValue(double x){
    HermiteFormula::get_DiffTable();
    double p = diffTable[0].get_value();
    double prod = 1.0;
    for (int i = 1; i < 2*PointNum; i++) {
        prod *= (x - inputPoint[(i-1)/2]);
        p += diffTable[i].get_value() * prod;
    }
    return p;
}